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We introduce a new method for calculating density perturbations in hybrid in- 
flation which avoids treating the fluctuations of the "waterfall" field as if they were 
small perturbations about a classical trajectory. We quantize only the waterfall field, 
treating it as a free quantum field with a time-dependent m 2 , which evolves from 
positive values to tachyonic values. Although this potential has no minimum, we 
think it captures the important dynamics that occurs as m 2 goes through zero, at 
which time a large spike in the density perturbations is generated. We assume that 
the time-delay formalism provides an accurate approximation to the density pertur- 
bations, and proceed to calculate the power spectrum of the time delay fluctuations. 
While the evolution of the field is linear, the time delay is a nonlinear function to 
which all modes contribute. Using the Gaussian probability distribution of the mode 
amplitudes, we express the time-delay power spectrum as an integral which can be 
carried out numerically. We use this method to calculate numerically the spectrum 
of density perturbations created in hybrid inflation models for a wide range of pa- 
rameters. A characteristic of the spectrum is the appearance of a spike at small 
length scales, which can be used to relate the model parameters to observational 
data. It is conceivable that this spike could seed the formation of black holes that 
can evolve to become the supermassive black holes found at the centers of galaxies. 
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I. INTRODUCTION 

Inflation [T] remains the leading paradigm for the very early universe. It naturally solves 
the cosmological flatness and horizon problems and is consistent with high precision measure- 
ments of the cosmic microwave background radiation [21 E]. Numerous models of inflation 
have been proposed, each adding features to the predictions of a scale invariant spectrum 
derived from single-field slow-roll inflation. Their motivation can be either some particle 
physics ideas coming from the standard model |I] or supersymmetric theories [5j [6], the 
need to explain some observation such as glitches in the CMB or supermassive black holes 
in galactic centers, or simply the extension of a theorist's toolbox in anticipation of the next 
set of high precision data, such as the upcoming Planck satellite measurements. 

Hybrid inflation was first proposed by A. Linde [7] and the name was chosen because 
this class of models can be thought of as being a hybrid between chaotic inflation and 
inflation in a theory with spontaneous symmetry breaking. The simplest hybrid inflation 
model requires two fields that we will call the timer and waterfall fields. The timer field 
corresponds the usual slow rolling field and is responsible for the scale invariant spectrum 
of perturbations observed in the CMB. The waterfall field is confined to its origin by the 
interaction with the timer field, giving a large constant contribution to the potential, which 
is also the main contribution to the energy density and hence the Hubble parameter. The 
potential governing the waterfall field changes as the timer field evolves, and at some point 
the minimum of the potential turns into a local maximum, and the waterfall field rolls down 
its tachyonic potential to its new minimum, where inflation ends. A characteristic feature 
of the density perturbation spectrum of hybrid inflation is the appearance of a large spike 
generated at the time when the waterfall potential turns tachyonic. The spike is generically 
at small length scales, and can potentially seed primordial black holes [8]. Primordial black 
hole formation and evolution has been studied in the past P-fl2lJ. but whether these black 
holes grow to become the supermassive black holes currently found in galactic centers is an 
open and intriguing possibility that we will address in a future publication. 

Usual inflationary perturbation theory is based on the study of quantum fluctuations 
around a classical trajectory in field space. However, in a purely classical formulation the 
waterfall field of hybrid inflation would remain forever at the origin, even after the waterfall 
transition, due to symmetry. It is quantum fluctuations that destabilize it and lead to the 
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end of inflation, so in a sense the classical trajectory has a quantum origin. Numerous papers 
have used various analytical approaches or numerical simulations to overcome this difficulty 
and approximate the spectrum of density perturbations [SJ El [T3TI2"T] . 

The method we use here has evolved from the early work in Kristin Burgess' thesis p3] , in 
which she studied a free-field model of the waterfall field in one space dimension, focusing on 
the time delay of the scalar field as a measure of perturbations. As in the model considered 
here, the waterfall field was described by a Lagrangian with a time-dependent m 2 , caused 
by the interaction with the timer field, m 2 evolved from positive values at early times to 
negative (tachyonic) values at late times. Such models are unnatural, since the potential 
is not bounded from below, but they nonetheless appear to be useful toy models, since the 
dynamics that generate the spike in the fluctuation spectrum occur during the transition 
from positive to negative m 2 . The evolution in the bottomless potential is realistic enough 
to give a well-defined time delay. Burgess studied the evolution of the waterfall field by 
means of a numerical simulation on a spatial lattice, using 262,000 points, calculating the 
power spectrum of the time delay by Monte Carlo methods. The method was slow, but 
for one choice of parameters she accumulated 5000 runs, giving a very reliable graph of the 
time-delay spectrum for this model. This line of research was pursued further in the thesis 
of Nguyen Thanh Son [14J, who repeated Burgess' numerical simulations with a new code 
(with excellent agreement). More importantly, Son and one of us (AHG), with some crucial 
input from private communication with Larry Guth, developed a method to short-circuit the 
Monte Carlo calculation. Instead of determining the power spectrum by repeated random 
trials, it was possible to express the expectation value for the random trials as an explicit 
expression involving integrals over mode functions, which could be evaluated numerically. 
The speed and numerical precision were dramatically improved. While Son's work was still 
limited to one spatial dimension, the possibility of extending it to three spatial dimensions 
was now a very realistic goal. In this paper we extend the calculation of the time-delay power 
spectrum in free-field models of hybrid inflation to three spatial dimensions, calculating the 
spectrum for a wide range of model parameters. 

In Section II we define the free-field model for the timer and waterfall fields that we will 
use to calculate fluctuations. We set up the equations of motion, define the notation of the 
mode expansion, and discuss the behavior of the mode functions. We make contact with 
a class of supersymmetric models that support hybrid inflation in Section III, presenting 
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the form of their potential and the range of parameters that they allow. Section IV gives a 
brief summary of the time delay formalism, and presents an approximation for calculating 
perturbations, developed earlier by Randall, Soljacic, and Guth. In Section V we develop a 
new method for calculating density perturbations in hybrid inflation that avoids any need to 
consider small fluctuations about a classical solution. Instead we show how the time delay 
power spectrum can be calculated essentially exactly in the context of the free field theory 
description. The result is given in the form of an integral over the modes which makes use 
of their known Gaussian probability distribution. In section VI we present an extensive set 
of numerical results over the parameter space of our model, where we are able to isolate 
the main factors that influence the density perturbation spectrum. In the limit of a light 
timer field, all quantities of interest are determined by the product of the timer and waterfall 
masses. We examine the models discussed in Section II as examples of realistic versions of 
hybrid inflation, and provide graphs showing the predictions of these models. Concluding 
remarks and directions of future work follow in Section VII. 



II. MODEL 



A. Field set-up 



Our first assumption is related to the expansion rate. We consider the metric to be 
exactly De-Sitter, even though this is only approximately correct. However, it changes only 
weakly during the slow roll inflation era, and we will terminate our calculation once the 
approximation loses its validity. Defining the Hubble constant during inflation as H, the 
scale factor is written as 

a{t) = e Ht (1) 
The model consists of two scalar fields. The "waterfall" field with lagrangian 



e 3Ht 



' 2 -e- 2 "'|V^-m^)|<^ (2) 



The usual 1/2 factors can be restored, if one writes = ^ (0i + ifa) where 0i and 02 are 
real scalar fields. The waterfall field must be complex, otherwise it will create domain walls 
as it rolls down from its initial value. The time-dependent mass of the field is controlled 
by a real scalar field, subsequently called the "timing" field. The important property of the 
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squared mass of is that it has to be positive initially and as ip evolves become negative. 
A general form is the following 

»*"=-4-(fy] 

We will choose r = 4 for most of our simulations. The lagrangian of the timing field is 



(3) 



e 3Ht 



T T v ' r/ 2""'' ^ 
The Lagrangians define the system up to an additive constant V in the potential, which 
is taken to be large enough, so that the variations in H are negligible during the era of 
interest. We neglected the interaction term from the Lagrangian of the timing field. This 
means that there is no back-reaction from the waterfall to the timing field. Physically this is 
a reasonable approximation before the waterfall transition, as well as afterwards, for as long 
as the waterfall field remains close to the origin. Mathematically, neglecting this term makes 
the equation of motion for the timing field de-coupled and in our quadratic approximation 
analytically solvable. 

Furthermore we do not examine perturbations arising from quantum fluctuations of the 
timing field. Before the waterfall transition they will give the nearly scale invariant spectrum 
that can be matched to the CMB observations. Apart from making sure that the long 
wavelength tail of the waterfall field perturbations does not contradict WMAP data, we will 
not consider these scales. After the waterfall transition the timer field perturbations will 
continue to be of the order of 10" 5 , hence they will be subdominant to the perturbations of 
the waterfall field by a few orders of magnitude, as we will see. The equations of motion are 

4> + 3H<p-e- 2Ht V 2 (p= -ml(t)(j) (5) 
i) + ZHj) - e~ 2Ht V 2 ^ = -m\i\) (6) 

If we take the timing field to be spatially homogenous, we get 

m = ^\ P = H^-l±^ 9 -- r ^ ] j ( 7 ) 

The value of the constant of integration was chosen so that ip(t) = tp c and rn 2 ^{t) = at 
t — 0. Both roots are negative, but the long time behavior is dominated by the larger of the 
two roots, which is 

'--"iVi/Tsii <*> 
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We will always choose -S- < 3/2 and not consider the case of a complex root. In fact, hybrid 
inflation models usually require the mass of the timing field to be well below the Hubble 
parameter, as in [5] and [5J. We choose to measure time in number of e-folds, hence we use 
iV = Ht. We rescale the masses similarly as fi^ = m^/H and \l§ = tuq/H. Furthermore the 
finite box size that we will use in our simulations is measured in units of jj and the field 
magnitude in units of H . We also define 

For a light timer field the reduced mass jl^ is proportional to the actual timer mass, jl^ = 

V 3 V H ) V 3^V>- 



B. Fast Transition 

Let's consider the speed of the transition. The transition happens at m| = 0. In order 
to quantify the speed of the transition, we will use the basic scale of our system, the Hubble 
scale. We will consider the transition duration to be the period for which \m^\ < H, meaning 
that the mass term in the equation of motion of the waterfall field is negligible. Assuming 
that fi^ > 1 we get 

±1 ^ (l _ e -^)^_l log g±l) (10) 

In the limit of Ji^ 3> 1 

AN (aw) 2 ^ 

Another measure of the speed of the transition is given by derivative of the waterfall field 
mass at iV = 0. 

1 dml(N) {n=o = M , ^ ^ ^ l_ (i2) 



R 2 dN « g^^y 

This shows that as long as the product fi^fi^p is somewhat larger than unity, the duration of 
the transition will be less than a Hubble time, meaning that the transition is fast! 



C. Mode expansion 

For purposes of our numerical calculations, we think of the universe as a finite box with 
periodic boundary conditions and a discrete spatial lattice. We choose the lattice to be cubic 
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with length b and Q points. This means that 

b 27T 

x = — I , k = —n , (13) 

where / is a triplet of integers between and Q — 1 and n is a triplet of integers between —Q/2 
and (Q/2) — 1. We can move between the finite discrete set of points and the continuous 
limit using the usual substitutions 

b\ 3 ^ f /2tt x 3 



£ £ , £ (14) 



QJ ^ J \t> 

x 

Our convention for the Fourier transform is 



k 



f(x) = Jd 3 ke fk -*f(k)=( 2 pj 

fc 

/(£) = (i)V d3x ^m = (^) 3 (^) 3 E e ^( f ) ■ ( 15 ) 

We will expand the waterfall field in modes in momentum space, 

i / n \ 3/2 

= (2^ W X^e^u&f) + S(k)e- ih3 u*(k,t)} , (16) 
which with Eq. ^ gives 

u(k, N) + 3u{k, N) + e' 2N k 2 u(k, N) = - e"^)^, iV) , (17) 
where k = jj- and an overdot denotes a derivative with respect to the time variable N = Ht. 

D. Solution of the mode function 

1. Early time behavior 

At asymptotically early times the k term dominates over the mass term provided that 
fi\ < 2. For r = 2 this is the case for fi^ < a/2, and for r = 4 it holds for m < a/5/4. These 
inequalities will hold throughout the parameter space of the models that we will examine, 
so we can neglect the mass term for N —¥ — oo. We then define a new function, following 
[22] as 

u(k,N)= 1 -J^e- 3N / 2 Z(z) , z = ~ke~ N . (18) 
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Neglecting the mass term in Eq. (17), we find 

2 d 2 Z dZ 
dz 2 dz 
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(19) 



which is the equation for a Bessel function of order 3/2. At early times the solution should 
look like a harmonic oscillator in its ground state, or equivalently the ground state of a 
free field in flat space, which is composed of negative frequency complex exponentials. This 
choice of initial conditions is the well known Bunch-Davies vacuum. For a review of scalar 
field quantization in de Sitter space and the corresponding vacuum choice see for example 
Ref. [23]. At early times the solution is given by 



u ~ 



(20) 



(21) 



where 

is a Hankie function, a linear combination of Bessel functions. (The phase is arbitrary, 
and the normalization is fixed by insisting that the field and the creation and annihilation 
operators obey their standard commutation relations.) Rewriting the original mode equation 
in terms of the new variable z, it simplifies to 

d 2 u 



2du 



dz 2 z dz 



+ u 



t4 



i - 



u 



(22) 



The k = mode is not captured by the procedure described here and is presented in detail 
in Appendix [A} 



2. General Solution 

We will now examine the general solution in a form that will be more appropriate for the 
numerical calculations that we have to perform. We can write the solution as 

u(k, t) = —l=R(k, t)e l0{lt) (23) 
V2kH 

and the differential equation separates in real and imaginary parts 

R - R6 2 + 3R + e- 2N k 2 R = ^(1 - e~^ N )R (24) 
2R8 + R6 + 3R6 = (25) 
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Integrating the second equation gives 



e -3N 



6 = const-^- (26) 



By comparing this with the early time behavior of the analytic solution 



u 1 e - N e fke - N (27) 



2HVk 

the phase equation becomes 

l p -3N 

e - -V < 28 > 

while the initial condition for the amplitude is given by the same asymptotic term to be 

R ->• e~ N (29) 
Inserting this expression in the equation for the amplitude function R 

~ k 2 e ~6N 



R - + 3R + e~ 2N ~k 2 R = - e~n N )R 

R A v 



(30) 



3. A closer look at the mode behavior 



Let us rewrite the equation of motion (Eq. 17) in a way that makes the time dependence 
of the solution more transparent 

Uk(t) + M k (t)+fi ff = , n 2 eff (k) = k 2 e- 2N +f4e~^ N -ft (31) 

We can distinguish different time windows with different behavior of the mode functions, 
based on the effective waterfall field mass. We will list these time windows here and then 
proceed to examine them one by one. 

1. N <^ 0, many efolds before the waterfall transition, in the asymptotic past 

2. N dev {k) < N < 0, a few efolds before the transition, where N dev {k) is the time at 
which a mode starts deviating significantly from the e~ N behavior, in particular starts 
decaying faster. 

3. < N < N tr (k) a few efolds after the transition, where N tr (k) is the time at which 
each mode starts growing. 
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4. iV 3> 0, the asymptotic future 

Now let us look at each of those time scales more closely. The asymptotic past is well 
described in the previous section and we see that all modes decay like e~ N . More precisely 
their magnitude behaves as \uk\ ~ \fik e ~ N ■ The first time scale Nd ev appears only for low 
wavenumbers. For N < we can keep only two of the three terms in the effective mass. 
Since /^e~^ 7V > fi 2 we will drop the pfy term, leaving the effective mass as ii 2 ^{k) = 
k 2 e~ 2N + jj^e~^ N . The time at which the two dominant terms become equal is 



Nde»(k) = -—z r ]og[- (32) 




2-^ 

For k > 11$ this time is not negative, hence we cannot drop \sh and our analysis fails. This 
transition, which happens only for k < fi^ signals a deviation of the behavior of the modes, 
which do not evolve as e~ N , but instead decay faster. 

Next we move to the actual waterfall transition time for each mode, which happens when 
the effective squared mass changes sign and becomes negative, or 

~ k 2 e~ 2N = nl(l- e~*l N ) (33) 

We will approximate the right hand side of the above equation with a piecewise linear 
function as follows 

4.-^")4f :(<1 « (34) 
For k < /i^e 1 ^^ the solution is found on the first branch and is 

Ntr(k) = \W (^] (35) 



where W is known as the Product Logarithm, or Lambert W function and is defined as the 
solution to the equation z = W(z)e w ^ z \ For small values of the wavenumber we can write 
the solutions as a Taylor series in k 

k 2 

N tr (k < v^^) = + Oik 4 ) (36) 

For k > /i^e 1 ^ we operate on the second branch and the transition time for each mode is 

N tr (k)= log (-) (37) 
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The behavior of the modes after the transition is different for different ranges of the 
timer field mass. If we consider the late time behavior of the mode equation, we can see two 
timescales introduced by the time dependent exponential coefficients. One is 0(1) and the 
other 0{l/jx^). We distinguish two cases: They can both be 0(1) or the second one can 
be larger than the first. The first timescale defines the time at which the equation becomes 
k- independent, meaning that all modes behave (grow) in the same way. The second time 
scale defines the time, after which the equation becomes time independent, meaning that 
after that all modes behave as pure exponentials. 

Let us first deal with the case of jl^ -C 1 meaning that the second time scale is much 
larger than the first one. Between the two timescales, that is 1 < N < the equation 

is independent of k 

R+3R = ^{l-e-~^ N )R (38) 

Since the evolution of the exponential term on the right hand side is by far slowest than the 
other timescales in the problem, we will treat 1 — e~^ N adiabatically. This leads immediately 
to the solution 



-3 + J9 + 4/4 (1 - e~~^ N ) 
R = R 0t MW , A(iV) = ¥ ^ (39) 

After a long time, this would mathematically settle to 



-3 + J9 + 4/xJ 
Ao = \ • (40) 

However, this is far beyond the time when inflation will have ended, hence it would physically 
never have time to happen (plus it is well outside the validity of our constructed potential). 

Let us choose = 10 and jfy = 1/10 to demonstrate our analysis. Some characteristic 
mode functions are presented in Fig. [T] 

We see that all modes behave similarly at late times, independent of their wave-number, 
as they should based on our late time analysis. Specifically, we can plot the ratio of the time 
derivative of each mode to its magnitude, as in Fig. |lj We call this the growth rate A = (J^J • 
We can see both phenomena. First, after N « 6 the modes behave identically. Second, the 
behavior of the mode approaches that of an exponential function (whose logarithm is a 
constant), but at a slower rate. In this example time needs to go on for several hundreds of 
efolds for the growth rate to set to a constant, which is calculated to be \(t — > oo) = 8.6119 
for /i^ = 10 and /x^ = 1/10. 
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FIG. 1: Mode functions for different comoving wavenumbers as a function of time in efolds. The 
model parameters are fj,^ = A and fi^ = 10. We can see the modes following our analytic 
approximation for the growth rate. Our analysis gives iV^e,, (1/256) « —7.9 and N tr (25Q) « 4.76, 
which are very close to the values that can be read off the graph. 

It is important to test our analytical approach to the late time behavior of the growth 
rate of mode function. As seen in Fig. [TJonce the mode functions evolve in a k - independent 
way, our simple analytical estimate for their growth rate is accurate to within a few percent, 
which gives us a very accurate expression for the growth rate and leads to the terms evolving 
as«~ e A (*)* ; where the time-dependent growth rate A(t) is slowly changing. 

As a test of our analysis, we can calculate the two important transition times 
N dev (1/256) = -7.9005 and iV ir (256) = 4.76457. We see that the calculated values agree 
very well with the behavior of the plotted modes. 

Let us briefly examine the situation where jl\ < 2. In this case the mode equation 
becomes k-independent and time independent at about the same time, that is a few e-folds 
after the waterfall transition. We choose = \ and /i^ = 1 and plot the results in Fig. 
[2j It is clear that the modes become both /c-independent and pure exponential (having a 
constant growth rate) at about the same time (N ~ 10). The asymptotic growth rate in 
this case is \(t — > oo) = 0.3028. 

Again we can calculate the two important transition times Nd ev { 1/256) = —7.8377 and 
iV tr .(256) = 5.39554, which agree once more with the behavior of the plotted modes. 
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FIG. 2: Mode functions for different comoving wavenumbers as a function of time in efolds. The 
model parameters are /i^, = k an d = 1. The horizontal line corresponds to the asymptotic value 
of the growth factor A. We can see how the mode functions reach their asymptotic behavior after 
10 efolds. Our analysis gives N^ev (1/256) —7.84 and A^ r (256) « 5.4, which are very close to the 
values that can be read off the graph. 

III. SUPERNATURAL INFLATION MODELS 

It is interesting to make contact between our abstract model and specific potentials 
inspired form particle theory. In general inflation models require small parameters in order 
to ensure slow roll inflation and produce the correct magnitude of density perturbations. It 
was shown in [5J, E] that supersymmetric theories with weak scale supersymmetry breaking 
can give models where such small parameters emerge "naturally" as ratios of masses already 
in the theory. We will not go into the details of such theories, but instead give the forms of 
the constructed potentials and use them as an application of our formalism. 

^MW(^/) + ^ + ^^ (4!) 

for what we will call model 1 and will be the primary focus of this work and 

V = M 4 cos 2 (c/>/V2f) + + X 2 — (42) 

2 4 

which we will call model 2. 

The first model can be taken with M' at one of three regions: the Planck scale, the GUT 
scale or an intermediate scale (~ 10 10 GeV). At each scale the rest of the parameters are 
adjusted accordingly to produce sufficient inflation and agree with CMB data. 
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Quadratic Approximation 


V 


m 


r 


tpc 




SUSY Model 1 


M 4 


M 2 
V2f 


4 


MV2M> 

Vf 




SUSY Model 2 


M 4 


M 2 
V2f 


2 


V2M 2 
/A 





TABLE I: Parameters of SUSY models and their counterparts in our quadratic approximation 



We will approximate the potential with a pure quadratic one with a time varying waterfall 
mass, of the form 



V(cj>^) = V -m 2 



2 ,1.2 



(43) 



where r = 4 for model 1 and r = 2 for model 2, as can be easily seen by the form of the 

interaction terms in both cases. The correspondence between the exact SUSY potential and 

our quadratic counterpart is shown in Table [Tj 

The parameters of the two models are restricted to fit CMB data, as shown in Fig. [3] 
To lowest order, in this potential dominated model, the Hubble parameter is constant 

and equal to 

ff=1 /5^ = J5^ (44) 



3 M n 



3 M„ 



End of Inflation 



In our simplified quadratic model inflation will never end. The waterfall field will roll 
forever down its tachyonic potential. However, we shall not forget that this is a mere 
Taylor expansion of more realistic potentials, which have a well defined minimum. We will 
use the supersymmetric potentials of E] as a concrete example to connect our purely 
quadratic potential to ones with more realistic shapes. In these supersymmetric models the 
potential has a cosine-like form and the minimum occurs at -J=-^ = |, where the inflaton will 
oscillate, terminating inflation and giving rise to (p)reheating. By making contact between 
the parameters of our potential and the physical parameters of the actual supersymmetric 
models, we can estimate the field value at which inflation ends. 

There are two strategies for defining </> e nd> the field value at which inflation ends. We 
can either pretend that the quadratic potential can be followed up to the end field value 
of the corresponding SUSY potential, or we can choose to end our calculation when the 
quadratic potential departs significantly from the actual SUSY potential that we are trying 
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FIG. 3: Parameter space for the two supernatural inflation models. The bottom right corner shows 
the parameter for model 2, while the other three show parameters for model 1, for different ranges 
of the mass scale M' 



to approximate. 

In the first case the end field value is at cn d = fn/y/2. To calculate the end field value 
for the latter case we will note that the cosine potential is accurately approximated by a 
quadratic as long as 0// I. We will call this ratio e and in this case we will end our 
calculations when e ceases being small. We can write these two cases in a unified manner, 
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as 

0cnd = ef (45) 

where e = tt/\/2 if we follow the quadratic potential all the way to the field value corre- 
sponding to the minimum of the SUSY potential and e < 1 if we stop our calculation at the 
point where the quadratic potential deviates significantly from the supersymmetric one. 

Using the values of the parameters taken from the supersymmetric models, we can esti- 
mate the end field value to be 

end ~ e 10 15 H (46) 

within one or two orders of magnitude for all cases of models considered in J5J E] ■ 

We will be using field values of this order of magnitude in our numerical calculations, 
whether we are dealing with the supersymmetric potentials or not. We will however ex- 
amine the effects of changing the end value of the field and show that it is minimal, easily 
understandable, and calculable. 



IV. PERTURBATION THEORY BASICS 



A. Time delay formalism 

The time delay formalism provides an intuitive and straightforward way to calculate 
primordial perturbations. Its basic principle is that inflation ends at different places in 
time at different times, due to quantum fluctuations. This leads some of the regions of 
the universe to have inflated more than others, creating a difference in their densities. The 
time-delay formalism was first introduced by Hawking [21] and by Guth and Pi [25], and 
has recently been reviewed in Ref. [26J. 

We will briefly describe the method here for the case of a single real scalar field. The 
universe is assumed to be described by a de-Sitter space-time, since the Hubble parameter 
is taken to be a constant. The equation of motion for the scalar field 4>(x,t) is 

* + 3H * = -^ + ^* < 47 > 

where the last term is suppressed by an exponentially growing quantity, so at late times it 
becomes negligible. We will omit the last term from now on. 
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We call the homogenous (classical) solution (f>o(t) and write the full solution, including a 
space dependent small perturbation 5(j) <C O as 

<f>(x,t) = 0oOO +5<f){x,t) (48) 

Plugging this into the equation of motion and working to linear order in 5(f) one can show 
that the quantity 5<j) obeys the same differential equation as <pQ. Furthermore the presence of 
a damping term implies that any two solutions approach a time independent ratio at large 
times. Thus, at large times we have (to first order in 8t) 

8<j>(x, t) -> -5r(f)0 o (t) => 0(f, t) -> <j) (t - 8r{x)) (49) 

This is the formulation of the intuitive picture of the time delay method. 

B. Randall-Soljacic-Guth approximation 

The usual calculation of density perturbations in inflation considers small quantum fluc- 
tuations around a classical field trajectory. In the case of hybrid inflation such a classical 
trajectory does not exist, since classically the field would stay forever on the top of the 
inverted potential. It is the quantum fluctuations that push the field away from this point 
of unstable equilibrium. One way to overcome this difficulty is to consider the RMS value 
of the field as the classical trajectory. This was done for example in [5] and [6] and is a 
recurring approximation in the study of hybrid inflation. 

Using the Bunch-Davies vacuum in the definition of the RMS value of the waterfall field 



rms = a/ (0\(J)(x, t)(J)*(x, t)\0) it it straightforward to calculate it using the mode expansion 

= ^~*' yx + d U<M, u l + = p E M*)i 2 ( 5 °) 

k , k> k 



The mean fluctuations are measured by 



A<f){k) 



■2ty 



<?xe ik - s (0(x)0*(O)) 



1/2 



2tt 



\ u k\ 



1/2 



(51) 



resulting in what will be called the RSG approximation for the time delay field 

^) *TiE«5(*)«iE(*) 
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FIG. 4: The time delay field calculated using the RSG formalism. The end time was taken to be 
15 e-folds after the waterfall transition and Li^p = ^ for all curves, while we varied ii^. 

There is an important comment to be made about the quantum mechanical nature of these 
density perturbations. In regular models of inflation quantum perturbations are scaled by 
h. We can think of them as modes with initial conditions that are of the order of h. The 
classical trajectory on the other hand does not have any quantum mechanical origin, hence 
does not scale with h. This means that in the limit of h — > the perturbations vanish, as 
one would expect will happen if one could "switch off" quantum mechanical effects. 

In the case of hybrid inflation on the other hand, what we call the classical trajectory 
(be it the RMS value or something else) is comprised of modes that originated as quantum 
fluctuations, hence is scaled by h itself. This means that even in the limit of h — > 0, the 
density perturbations in hybrid inflation remain finite! By explicitly restoring h in the 
formulas of the paper, the reader can formally arrive to the same conclusion. 

Some plots of the time delay field calculated using the RSG approximation are shown in 
Fig. |4j The reduced mass of the timer field was taken to be = ^ while we varied the 
waterfall field mass. We fixed the time at which inflation ended to be 15 e-folds after the 
waterfall transition. 

V. CALCULATION OF THE TIME DELAY POWER SPECTRUM 

The usual method to calculate the primordial perturbation spectrum would involve either 
making some approximations (more or less similar to the RSG) or using a Monte Carlo 
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simulation. The former suffers from the lack of a classical trajectory that invalidates the 
usual perturbation method, while the latter would be computationally costly in three spatial 
dimensions. We will therefore implement an alternate method that reduces the calculation 
of the spectrum of the time delay field to the evaluation of a two dimensional integral and 
does not need a classical trajectory to do so. 



As discussed at the end of Sec. IIP 3 the behavior of the mode functions at asymptotically 
late times (t —> oo) is given by 



u(k,t oo) ~ e Xot u(k) , (53) 



where Aq is given by Eq. (40). If we define for all times 



• V R(k,t)R(k,t) 

X(t) = = m , (54) 



'k 2\k\ 

then at late times \{t) — > Ao- Since \(t) changes very slowly, we can take it as a constant 
around the time of interest. 

To discuss fluctuations in the time at which inflation ends, we begin by defining to as the 
time when the rms field reaches the value e nd ; which we have chosen to define the nominal 
end of inflation: 

0rmsOo) = 0end ■ (55) 

Since at late times all modes, to a good approximation, grow at the same exponential rate 
X(t), we can express the field t) at time t = t Q + St in terms of the field 4>(x, to) by 

|0(f,t)| 2 = |0(x,t o )|V A5 *. (56) 

If t is chosen to be the time tend(^) at which inflation ends at each point in space, then 
(f)(x, t en d(x)) = 4>end = 0rm S (^o), and the above equation becomes 

<f>L s (to) = \mto)\ 2 e 2X5 \ (57) 

which can be solved for the time delay field 8t(x) = t end (x) — t : 

St{x) =2xHl^M) ' (58) 



Rescaling by the rms field 



(59) 
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we can write 



St(x) = ^log|0(f,t o )| 2 



Using this expression, we can write the two-point function of the time delay field as 

5t{x)6t(0)) = ^ (log t ) | 2 log |0(O, t ) I 2 



(60) 



(61) 



which can be evaluated, since the probability distributions are known. To continue, we can 
decompose the complex scalar field in terms of the real fields Xf 



{x, t)=X 1 + iX 2 , 0(0, t) = X 3 + iX A . 



(62) 



The average value of a function F of a random variable X with probability distribution 
function p(X) is given by 

(F[X}} = [ dXp(X)F[X] . (63) 



Since this is a free field theory, we can take the four random variables Xi(x) to follow a joint 
Gaussian distribution with 



p(X) 



1 



, exp l --X'H-'X 

(27r) 2 v/det(S) V 2 



1 



V,-- 



(XiXi) 



(64) 



A function of the X[s then has the expected value 

A 

1 



(F[X)} 



JJdXi 



1 



(2^Vd^) eXP l-I XE X)FlX] 



(65) 



The new fields Xi can be written in terms of the original complex field (p as 

1 



X 1 

x 3 



(j)(x) + (j)*(x) 
0(O) + 0*(O)" 









4>{x) — <p*(x) 


2i 





x d 



2i 



0(O)-0*(O) 



(66) 



The components of the variance matrix E can be easily calculated using the commutation 



relations for the creation and annihilation operators in cj)(x, t), from Eq. (16). Due to the 
high degree of symmetry the matrix itself has a very simple structure: 



( \ A \ 

\ A 

A \ 

V0 A \) 



(67) 
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where 



A(x,t ) = (X 1 X 3 ) = (X 2 X 4 ) = l - (0*(f,t o )0(O,t o )) = Wk,t )\ 2 e rk -* , (68) 



where 



u(k, i) = 



u(k, t) 



(69) 



0rms(i) 

Since u(k,t) actually depends only on the magnitude of the wavenumber, because of the 
isotropy of the problem, we can do the angular calculations explicitly in A and leave only 
the radial integral to be calculated numerically. Then 

1 1 



St(x)St(0) 



4A 2 (27r) 2 [i-A 2 ] 



x exp 



dX ± dX 2 dX 3 dX 4 log(X 2 + X 2 ) log(X 2 + X{) 



4[| - A 2 ] 



X\ + X\ + X\ + X\ - 4(XiX3 + X 2 X 4 ) A 



(70) 



Changing to polar coordinates 



X\ = T\ cos 6\ , X 2 — r\ sin 6*i 
X 3 = r 2 cos 6 2 , X 4 = r 2 sin 6 2 , 



(71) 



the integral becomes 



ttA 2 (1 -4A 2 

r i + r i ~~ 4Arir 2 cos 6* 



;>27r ;>oo i>oo 

i d6 / rxrfri / r 2 dr 2 log(ri) log(r 2 ) 
o Jo Jo 



x exp 



1 -4A 5 



(72) 



where we redefined the angular variables as 9 = 9\ — 6 2 and 9 = 9i + 9 2 and integrated over 
9. Changing also the radial variables 



ri = r cos <p , r 2 — r sin i 



(73) 



5t(x)5t(6) 



^ t>2-K poo 

— — — — / d9 / sin 20 / drr 3 log(r cos0) log(r sin0) 

7TA 2 (1 - 4A 2 ) J J Jo 



x exp 



(1 -2Asin20 cos6>) r 2 " 
1-4A 2 _ 

The radial integration can be performed analytically 



dr r 3 log(ar) log(6r)e cr = 



(74) 



(75) 



o 
1 

8? 



7T 



( 7 - 2) 7 + — - 21og(a6)( 7 - 1 + log(c)) + 41og(a) log(6) + log(c)(2 7 - 2 + log(c)) 
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where a = cos0, b = sin0, c = (jz^M) 0- ~ ^A si 11 20 cos 8) and 7 is the Euler constant 
7 « 0.57721. 

Finally, the spectrum of the time delay field is defined by 



Sr(k) 



k^ 3 



( 1/2 

d 3 xe ik -* (8t(x)8t(0) 



(76) 



2tt 

Calculation in the two limiting cases x — > and x — > 00 (or x — >■ 6 in our case) can be 
done analytically. 

1. For x — y several terms in the integral diverge, since A — > |. In this case we have 
only two degrees of freedom instead of four, since we consider a complex scalar field 
at one point in space. The integral becomes 

[St(0)6t(6)) = ±;J d -^^e-^ + ^\o g \Xl + Xl) = 

/2n poo 1 / 2 \ 

d6 ^ dr re~ r2 log 2 (r 2 ) = — ^ 7 2 + - J . (77) 



4ttA 2 

2. The x —7- 00 limit is much easier to handle. We recognize that A(x) is simply the 
Fourier transform of \uk\ 2 - Since is smooth, A(x — > 00) — >■ 0, and therefore 5t(cxo) 



is uncorrelated with <5t (0) . Eq. (70) can be seen to factorize, giving (St (00) 8t(Q) 
(^5t(6)y, where 



1 

~7rA 

7 
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5t(0)) = ^ y rfX^Xa log(X 2 + X 2 ) exp [-(X 2 + X 2 2 

/■°° _ 2 

d9 I r dr log re r 
Jo 

(78) 



2A . 

Combining these results, we see that the probability distribution for 5t(0) has a standard 
deviation a — \ j (5t(0) 2 \ — (st(6)\ = 7r/(2y / 6A). While the first limit above is needed for 



programming the numerical calculations, since the integral of Eq. (70 ) cannot be numerically 
evaluated at x — 0, the second limit can be used as a numerical check. 

The same method can be applied to the exact calculation of any higher order correlation 
functions. Especially the non-Gaussian part of the power spectrum f^L can be read off 
from the momentum space Fourier transform of the three-point correlation function in posi- 
tion space {5t( y Xi)6t(x 2 )St(x 3 )) = (5t(xi)St(x 2 )St(0)) . Taking the Fourier transform we can 
compute (5t(ki)5t(k 2 )5t(k 3 )\, from which we can extract the properties of the bispectrum. 
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(tft(gi)ft(g 2 )ft(0)) = y / ^.,i//-2 / sin0d0 / #F( 7l , 7 2 ,0,0) (79) 



The form of the three point function in position space is 

\2 fl-K fTT r 2~ 

A 3 " 

where ^(71,72, 0, 4>) is a function of four angular variables. Calculations regarding the form 
of the bispectrum will be published elsewhere. 

VI. NUMERICAL RESULTS AND DISCUSSION 

Let us begin by plotting one example of the free field theory (FFT) calculation of the 



time delay power spectrum, from Eqs. (74) and (76), along with the corresponding curve 



derived using the RSG approximation, Eq. (52). We use the sample parameters = ^ 



and fie/, = 20. Both calculations give a spike, but a spike of different width, different height 
and different position. Let us rescale the RSG result as follows 

<57RSG,rescaledO) = A Strsq (Bk) , (80) 

where A and B are 0(1) constants calculated by requiring the peaks of the FFT and RSG 
distributions to match in position and amplitude. The results are plotted in Fig. |5j We 
can see that the FFT and RSG curves do not seem similar. However the rescaled RSG 
curve seems to follow the FFT curve very well, as was first noticed by Burgess [13J. Based 
on our simulations the curves generally tend to agree better for low wavenumbers, up to 
and including the peak, and start deviating after the peak. The rescaling parameters vary 
with the field masses chosen and for the particular choice of Fig. [5] were calculated to be 
A = 0.6152 and B = 3.25 . We do not yet fully understand this behavior, but we are 
studying it both analytically and numerically and will present our findings in a subsequent 
paper. 

We will now do an extensive scan of parameter space {/^, /A/>} in order to have reliable 
estimates on the magnitude and wavelength of the perturbations. This is important both 
to make sure that CMB constraints can be satisfied as well as to study the formation of 
primordial black holes that might lead to the supermassive black holes found in the centers 
of galaxies. Since the original motivation for this paper has been the supersymmetric models 
first presented in [5] and [B], we will present the results for the perturbations in these models. 
However, our quadratic approximation holds for more general hybrid inflation models. Hence 
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FIG. 5: Comparison between the RSG and FFT methods. The end time was taken to be 15 e-folds 
after the waterfall transition and = 20 and fj,^ = 1/20. We can see that the spectrum of the 
time delay field calculated in the free field theory agrees very well with the rescaled version of the 
RSG approximation A5trsg (Bk). 

it is important to make a model-independent parameter sweep. This will provide a more 
general set of predictions of this class of models. We will give both exact power spectra, as 
well as try to isolate the dominant features and provide a qualitative understanding of their 
dependence on the model's parameters. 

A. Model-Independent Parameter Sweep 

There are several model-dependent parameters that give us some control over the proper- 
ties of the resulting power spectrum. Initially we will fix the value of the field at the end of 
inflation to be |0 e nd| = 10 14 in units of the Hubble parameter. With this assumption (which 
will be relaxed later), we can calculate the properties of the power spectrum as a function 
of the masses. Initially we fix the reduced timer field mass to be = ^ and vary the mass 
of the waterfall field. The results are shown in Fig. [6] We have plotted (clockwise from the 
top left) 

1. The end time of inflation, defined as the time when the RMS value of the field reaches 
the end value. 
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FIG. 6: Parameter sweep for constant timer field mass /i^ = 1/20 and constant end field value 
0end = 10 . Data points are plotted along with a least square power law fit. The same trend 
is evident in all curves. The time delay spectrum grows in amplitude and width and is shifted 
towards larger momentum values as the mass product decreases. Also inflation takes longer to end 
for low mass product. 

2. The maximum amplitude of the spectrum of the time delay. 

3. The comoving wavenumber at which the aforementioned maximum value occurs. 
Thinking about black holes, this is the scale at which black holes will be most likely 
produced. 

4. The width of the time delay distribution in the logarithmic scale, taken as Ak = 
log 10 where k±i/ 2 are the wavenumbers at which the distribution reaches one 
half of its maximum value. 
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We see that all the plotted quantities show a decreasing behavior as one increases the mass 
product. In order to quantify this statement, we fitted each set of data points with a 
power law curve of the form y = ax b + c. The scaling exponent b for the various quantities 
was b tcnd ft -0.88, b 5rmax ft -0.34, b kmax ft -3.219, b Ak ft -1.17. As a comparison, the 
corresponding best fit exponent of the growth rate A as a function of the mass product is 
b x ft -0.85. 

In order to get a better understanding of what these parameters actually mean, we plot 
three characteristic spectra for three values of the mass ratio in Fig. [7j We also rescale the 
spectra by the growth factor A. This probes the actual form of the two point correlation 
function, as seen in momentum space. That is, it shows the evaluation of the spectrum in 



Eq. (76), while ignoring the factor of 1/A 2 in the evaluation of (5t(x)5t(6)j from Eq. (74). 



Before continuing to a more thorough examination of parameter space, let us understand 
how changing the field value at the end of inflation will change our results. Fixing the 
product of the reduced masses equal to 2 (//^ = ^ and //^ = 40), we let the field value 
</w vary by four orders of magnitude. The results are shown in Fig. |8j It is evident that 
the curves for Sr(k) are of identical form and slightly different magnitude. The last graph 
shows the product A ■ Sr(k) for the two curves at end = 10 12 and en( j = 10 16 , plotted 
respectively as a green thick and a black thin line. It is seen, that once rescaled the two 
curves fall exactly on top of each other, meaning that the actual integral that gives us the 
two point function in position space is time independent, once we enter the region where 
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FIG. 8: Perturbation spectrum for varying field value at i e nd for constant masses. The time delay 
curves are identical in shape and differ only in amplitude. This variation is entirely due to the 
different value of the time dependent growth factor A, which differs for each case because inflation 
simply takes longer to end for larger end field values. 

all modes behave identically. Furthermore if one takes the product of the maximum value 
of the time delay times the growth parameter {8r max ■ A) for the different values of e nd the 
result is constant for the range explored here to better than 1 part in 10 6 , meaning that 
they are identical within the margins of numerical error. Thus, changing the value of the 
field at which inflation ends can affect the resulting perturbation spectrum only by changing 
the growth parameter A, for which we have a very accurate analytical estimate in the form 
of Eq. (39). From this point onward, we will keep the end field value fixed at e nd = 10 14 



and keep in mind that the fluctuation magnitude can change by 10% or so if this field value 
changes. 
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Once we fix the field magnitude at the end of inflation we have two more parameters 
to vary, namely the two masses: the actual timer field mass and the asymptotic tachyonic 
waterfall field mass. The two masses can be varied either independently on a two dimensional 
plane or along some line on the plane, in a specific one-dimensional way. Fixing one of the 
two masses is such a way of dimensional reduction of the available parameter space, as we 
did before. Another way to eliminate one of the variables is to fix the mass product and 
change the mass ratio. This will prove and quantify the statement, that (at least for heavy 
waterfall and light timer fields) the result is controlled primarily by the mass product. 

We fix the mass product at = 2. The results are shown in Fig. [9] The curves are 

of identical form and everything is again controlled only by A. On the top left figure we 
plotted A St for the two extreme values and the curves fall identically on top of each other 
(color-coding is as before). Furthermore if we calculate the product A • ST max for different 
values of the mass ratio we get a constant result 0.1225 ± 2 • 1CT 5 where the discrepancy 
can be attributed to our finite numerical accuracy. The second feature of this calculation is 
the extremely flat part of the end-time, growth factor and maximum time delay curves for 
large values of the mass ratio and the abrupt change as the mass ratio gets smaller. For the 
value of the mass product that we have chosen, this transition happens as the timer field 
mass approaches unity. Let us look at the expansion of the effective waterfall field mass 

»leff = A*J (l - e~^ N ) = nffiN{l - %N + ...) (81) 
When the second term in the expansion cannot be neglected, the dynamics of the problem 
stops being defined by the mass product alone. This explains the abrupt change we see as 
we lower the mass ratio. By doing the same simulation for different values of the fixed mass 
product we get similar results. 

We can now do the opposite, that is fix the ratio and change the mass product. The 



results are shown as the open circles in the top two diagrams of Fig. 10, and in the lower 
diagrams of the figure. There are two main comments to be made. First of all, in the 
case of a fixed ratio, the growth rate A does not solely determine the results. Rescaling the 



spectrum by A not only fails to give a constant peak amplitude (Fig. 10, lower left), but the 



result are spectra of different shapes (Fig. 10, lower right). On the other hand, the data 
points taken with a constant mass ratio and a constant timer field mass (//^ = ^), as shown 



by the +'s on the upper diagrams of Fig. 10, fall precisely on the same curve! This clearly 
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FIG. 9: Fixing the mass product at 2 and varying the mass ratio. There is significant variation 
only for low mass ratio, when the light timer field approximation loses its validity. Furthermore the 
curves of maximum time delay amplitude and 1/A follow each other exactly up to our numerical 
accuracy. Finally by rescaling the spectra by the growth factor A they become identical for all 
values of the mass ratio. 

demonstrates that the only relevant parameter, at least for a light timer field, is the mass 
product! 

We see that contrary to the fixed product case, the results for fixed mass ratio do not 
depend solely on A. We have established that the the most important factor in determining 
the time delay field is the product of the waterfall and timer field masses, especially for a 
light timer field. 
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FIG. 10: Fixing the mass ratio at 900 (open circles) or the timer mass at = 1/20 (+'s). 
The time delay spectra for different mass products show no common shape characteristics and 
remain different even when rescaled by A. Furthermore the end time and maximum perturbation 
amplitude curves are identical for constant mass ratio and constant timer field mass, proving that 
indeed the mass product is the dominant parameter. 

B. Supernatural Inflation 

We now turn our attention to the supernatural inflation models that were studied in [5] 
and [6]. We will examine each of the four cases separately. 



Let us start with the first SUSY model (described by Eq. (41)) with the interaction- 
suppressing mass scale M' set at the Planck scale. The mass of the timer field was calculated 
to be 50 to 100 times less than the Hubble scale, while the asymptotic waterfall field mass 
was more than 20 times the Hubble scale. This means that the model is well into the 
region where the two masses are separated by a few orders of magnitude. According to the 
analysis of the previous section, we expect the mass product to be the dominant factor in 
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FIG. 11: First Supernatural inflation model with M' at the Planck scale. The spectra corresponding 
to the maximum and minimum mass product are shown. We observe good agreement with the 
results of the model independent parameter sweep of the previous section, because the timer field 
mass is much smaller than the Hubble scale. 



the generation of density perturbations. In the left part of Fig. [TT] we see the mass product 
for this model. We can see that the mass product varies less than 15%. It is hence enough 
to calculate the time delay spectra for the two extreme values and say that all other values 



of the mass product will fall between the two, as shown in 11 



Putting the mass scale M' of the first SUSY model at the GUT scale changes the masses 
as well as the Hubble scale by one order of magnitude. However the reduced masses and 



their product have very similar values as before. This is shown in Fig. 12 



It is worth noting that these two SUSY models contain a very light timer field, hence the 
results should be the same as our previous parameter space sweep with a constant light timer 



field. If one compares Fig. 11 and Fig. 12 with Fig. km we indeed see excellent agreement 



for the amplitude and width of the time delay spectrum. 

When setting the mass scale M' at some lower scale of 10 11 GeV, the reduced 
timer and waterfall masses become 0(1). This means that in this case the parame- 
ter A saturates faster and the perturbation spectrum reaches its asymptotic limit ear- 
lier and becomes time-independent from that point onward. Furthermore the actual 
value of the growth parameter A is smaller, leading to an enhanced perturbation am- 
plitude, the largest among the models studied here. The mass product changes by a 



factor of 2.5 as seen in Fig. [13} We choose five points in the allowed interval of mass 
values and calculate the corresponding curves. The specific values of the mass parame- 
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FIG. 12: First Supernatural inflation model with M' at the GUT scale. The spectra corresponding 
to the maximum and minimum mass product are shown. There is again good agreement with the 
results of the previous section. 
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FIG. 13: First Supernatural inflation model with M' at the intermediate scale. Five representative 
pairs of masses were chosen and the corresponding time delay curves are shown. This model can 
give maximum time delay of more than 0.1. 

ter M are M = 1.06 • 10 10 GeV, 2.4 • 10 10 GeV, 5.42 • 10 10 GeV, 1.23 • 10 11 GeV, 2.77 • 
10 11 GeV. The corresponding pairs of reduced waterfall and timer masses are {/i^,, 1//^} = 
{3.19,4.48}, {2.58,2.97}, {2.22,2.05}, {1.99,1.47}, {1.83,1.09}. The points on the mass 



product graph are color coded to match the corresponding time delay curve in Fig. 13 We 
can see that since the mass products have a larger variation, the resulting spectra have quite 
different time delay spectra. Also, since the timer is not much lighter than the Hubble scale, 
the curves do not scale according to our previous analysis. 



We finally consider SUSY model 2, Eq. (42), with the ijr<jr interaction term. Again the 
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FIG. 14: Second Supernatural inflation model. Three representative pairs of masses were chosen 
and the corresponding time delay curves are shown. 



reduced masses are 0(1), so we expect a small A leading to a large amplitude perturbation 
spectrum. The mass product varies around 1 by less than ±15%. We choose three values of 
the mass product (the two extrema and an intermediate one) and plot the resulting curves 



in Fig. 14 The specific values of the mass parameter M are M = 1.080 • 10 GeV, 1.006 • 
10 11 GeV, 9.376T0 11 GeV and the corresponding pairs of reduced waterfall and timer masses 
are {/v, 1//^} = {2.697, 2.367}, {2.330,2.262}, {2.007,2.170}. 



VII. CONCLUSIONS 



We presented a novel method for calculating the power spectrum of density fluctuations in 
hybrid inflation, one that does not suffer from the non-existence of a classical field trajectory. 
We used this method to numerically calculate the power spectrum for a wide range of 
parameters and concluded that in the case of a light timer field, all characteristics of the 
power spectrum are controlled by the product of the masses of the two fields. In particular 
the amplitude was fitted to a power law and found to behave as 5T max ~ O.O3(yU0yU,/,)~ ' 34 
and the width in log-space as Ak ~ 1.7(/i0 / u^,)~ 117 . Furthermore we made connection to 
SUSY inspired models of hybrid inflation and gave numerical results to their power spectra 
as well. For the SUSY models with a light timer field the numerical results were in excellent 
agreement with our fitted parameters. 

Work is currently under way in refining and extending the formalism. Understanding 
the rescaling properties between the RSG approximation and the exact result could provide 
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further insight into the physics of the problem and provide quasi-analytical approximation 
of well controlled accuracy. We will also apply our results to estimating the number and size 
of primordial black holes and try to make contact with astrophysical observations regarding 
supermassive black holes in galactic centers. Finally we are examining the predictions of 
our model for the non-Gaussian part of the perturbation spectrum. 
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Appendix A: Zero mode at early times 

— * 

The k — mode is not captured by the procedure described in the main text. If we 
consider this mode alone for asymptotically early times, so that we keep only the exponential 
in the mass term 

u + 3m = -/i|e" A ^M (Al) 

This can again be solved in terms of Bessel functions by defining a new variable and a new 
function as 

z = ae~ ii * N/2 ; u (Q,N)=z p Z(z) (A2) 

The mode function becomes 

?g +i f( 1 + 2fi -$) + z(p-W + 4g) =0 (A3 ) 

dz 2 dz y n\) y n\ a 2 ^ J 

The standard form of the differential equation that gives Bessel functions is 

f^ + z dZu + ( 2 ?_j )Zv = (A4) 
dz 1 dz 
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By appropriately choosing the constants a,/3 and v the two equations can be made identical. 
The choices are 

= 4> « = ^> "=4 (A5) 
^ip ^ip 

Finally introducing an arbitrary constant of normalization No, the solution for the zero mode 
at asymptotically early times becomes 



^,N) = N e^ N ' 2 H^\z) , z= 2 ^±e-^ N/2 (A6) 

Pip 



The normalization factor N can be defined using the Wronskian at early times. The Wron- 
skian at all times is defined as 

wCk, t ) = uCk,t)^^-^pu-(-k,t) (A7) 

Taking the time derivative and using the equation of motion 

= „(M)^zM - ^ U .(-M) = -,W(M) 

=>W(k,t) = f(k)e- 3Ht (A8) 

— * 

Since f(k) is by definition independent of time, we will compute it at approximately early 
times, where we know the solution in analytic form and the solution is 

W (k^0)= ie ~ 3N , W(k = 0) = 2 ^^N'e- 3N (A9) 

pi 

Requiring that the Wronskian be a continuous function of k at all times we can extract the 
value of N . 

^/3S (A10) 

Appendix B: Initial Conditions 

We can rewrite the mode equation as a system of three coupled first order differential 
equations. 

d0 _ ke-™ 

dN ~ — ffT" (B1) 
§ = ^ -3R- e-^-R + - e-*l»)R (B3) 



36 



In this notation, R is one of the three independent functions. 

This is not a system of three coupled ODE's in the strict sense. We can first solve the 
two equations j§ and ^ as they do not contain any terms involving 9 or its derivative. We 



can then integrate forward in time, using the calculated values of R(N). Furthermore 
it is clear that the equations only depend on the magnitude of the wavenumber, as was 
expected due to the isotropy of the problem, so we need only solve the mode equations for 
one positive semi axis. 

We know from the analytical solution at early times that 



R(N -oo) ->■ e 



—N 



(B4) 



Since we have to start the numerical integration at some finite negative time without losing 
much in terms of accuracy, we refine the initial condition by including extra terms in the 
above expression. We will then start numerically integrating when our expansion violates 
the desired accuracy bound. We define the correction to the asymptotic behavior as 5R(N) 
such that 

R(N) = e~ N + 5R(N) (B5) 
We will expand 5R in powers of /i^ and e N . 



5R = 



,3V 



,5iV 



2k 2 8k 4 lQk 6 

IN 



4k 2 



1 - e-^ N 



,3JV r 



+ 



16A; 4 



4 + e-^(/2j - 6/$ - 4)] 



64A; 6 



86 + e -W fa* _ 14*. + 53^ _ 25 ~2 _ g6) 



3 3V 



32k 4 



1 - e-^ 1 



64k 6 



29 - e" A * JV (9/iJ - 65//^ + 58) 



+e -2^JV( 14 ^4 _ 65 ~^,2 + 2g ^ j + 15fl 



5V 



128/fc 6 



1 - e'W 



(B6) 



We can now choose the initial expansion for R(N) to calculate the expansion for the phase 
0(N). 

The asymptotic behavior, given by the standard definition of the Hankel functions is 



6{N ->■ -oo) = ke~ N - 7T <=> 0(N ->■ -oo) + vr = fee 



(B7) 



We define = 0(iV) + 7r =>- # = in order to keep track of the constant phase factor without 
carrying it through the perturbation expansion. 
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FIG. 15: Mode functions for = 22 and mu^p = 1/18. The left column is calculated for k = 1/256 
and the right for k = 256 

Defining the corrections to the early time behavior of the phase as 

= e - N [k + 59(N)} (B8) 

we can construct a similar expansion as the one for SR. Although our formalism does not 
require knowledge of 0(N), we included it for completeness. 
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